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<N : 
^3 : 

■ We have used the vectorised and parallehsed magnetohydrodynamics code 
NIRVANA on the NEC SX-5 and the new SX-6 instaUation in parallel mode 
to simulate the interaction of jets with a galactic wind that might be typical 

I for the star-bursting radio-galaxies of the early universe. 
T-H ■ The two simulations, one axisymmetric and one in 3D show that the jet 

\ pierces and destroys a thin and dense shell produced by the pre-installed 

■ super wind. 

I We suggest that small radio galaxies at high redshift might be absorbed on 

QQ ■ the blue wing due to the galactic wind shell, and possibly on the red wing due 

I to a cooling flow. In larger sources the jet cocoon will fill the wind cavity and 

'nI" ■ accelerate the shell. The Ray leigh- Taylor instability will then disrupt the shell 

' and disperse dense, possibly star-forming fragments throughout the region. 

Q ■ 1 Introduction 

^: 

c/3 ■ Radio galaxies are associated with an enormous release of energy, exceeding 

; 10^° Watts (Ghisellini 2003; Krause & Camenzind 2003). From the power 

^ ■ house in the center of certain galaxies, where supermassive black holes of 

billions of solar masses are believed to reside, the power is channeled through 

■ narrow jets out to distances of several million light years, far outside of the host 
' galaxy. They stimulate emission throughout the electromagnetic spectrum 

that can be observed up to the largest distances of the visible universe. Being 
the likely progenitors of today's brightest cluster galaxies (Carilli et al. 2001), 
radio galaxies at high redshift mark the peaks of highest density in the early 
universe, and highlight the formation of galaxy clusters. 

At redshifts of z > 0.6 extended optical continuum and emission line re- 
gions, aligned with the radio structures, become prominent companions (Mc- 
Carthy 1993). Beyond a redshift of z w 2, most of the optical emission origi- 
nates from the Lyman a transition of hydrogen. Huge Ly a halos have been 
observed in the young universe (e.g. Reuland et al. 2003). Their size often 
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exceeds the radio size. Only for radio extents smaller than 50 kpc^ , the Ly a 
emission, that has turbulent velocities of typically 1000 km/s, is totally ab- 
sorbed with a velocity width of several 10 km/s (van Ojik et al. 1997). The 
absorbers are typically blue-shifted (van Ojik et al. 1997; De Breuck et al. 
2000). 

Optical emission lines are produced by clouds at temperatures of 10"^ to 
10^ K. These are thought to be embedded in thin and hot plasma in pressure 
equilibrium. The emitting gas is cither ionised by ultraviolet radiation, e.g. 
from hot stars or the central quasar, or coUisionally excited due to shocks 
associated with the jet expansion. The location of the absorbing screen is still 
a matter of debate. Based on the determination of differing mctalicitics in the 
emitting and absorbing gas components, some authors have suggested an ex- 
tended low density shell around these systems (Binette et al. 2000; Jarvis et al. 
2003; Wilman et al. 2003). This has the disadvantage that the low velocity 
width of the absorbers seems to be hard to explain. Based on hydrodynamic 
simulation, one of the authors has suggested a different model (Krause 2002). 
According to that, the emission line region is surrounded by a high density 
shell. This shell is formed when the gas behind the leading bow shock is dense 
enough to cool on the propagation time-scale. This part of the gas looses pres- 
sure support and is compressed in a thin and dense shell that would produce 
the observed absorption. The latter model may produce too much luminosity 
in the Ly a and X-ray channel. The combination of a jet with a galactic wind 
is similar to the that model with the exception that lower luminosities of the 
absorbing shells are involved. 

The jet-galactic wind interaction scenario is detailed in sect. 2, followed 
by a discussion of the computational aspects in sect. 3. We present the results 
of the simulations in sect. 4. They are discussed in sect. 5. 

2 Details of the interaction of an extragalactic jet with a 
galactic wind. 

This model (compare Fig. 1) starts with a starburst. This means that the stel- 
lar winds and supernovac inject gas and energy into the interstellar medium, 
which leads to an outflow of gas. The outflow is headed by a bow shock 
(Fig. la). After the cooling time, the material behind this shock will cool and 
form a dense shell (Fig. lb), which may absorb emission from gas inside of the 
shell. The shell may be Rayleigh- Taylor unstable and some gas may fall back. 
Additionally, cjecta from the stars of the galaxy fill the cavity. When the jet 
starts (Fig. Ic), the gas inside of the cavity may be excited to emit strong 
Ly a due to photoionisation by the central quasar, which has to be assumed 
to be active when it ejects jets, or by shocks associated with the jet expansion. 
This emission will be absorbed on the blue wing due to the expanding shell, 

^ kpc = 3260 light years = 3 x 10^^ m 
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1. Starburst causes galactic wind 



2. Thin shell forms after cooling 




3. Jel inside the Ihin and dense shell 



4. Jet dis rupts t he shell 

^SESIffiB 




Fig. 1. Sketch of the interaction of a jet with a galactic wind adopted here. First, the 
supernovae associated with a starburst produce a galactic wind (a) . After the cooling 
time, the gas inside of the bow shock cools and forms a dense shell (b). The shell 
may become Raylcigh- Taylor unstable and loose some material to the inside. Then, 
the jet starts (c). The stars still eject gas which forms an emission line region by 
photoionisation or jet induced shocks. The thin and dense shell absorbs at a narrow 
velocity range. Wo suggest the small radio galaxies to be in this phase. When the jet 
pierces the shell and the exhaust jot plasma fills the cavity, the shell is accelerated 
because of the high pressure in the interior (d). The shell is then disrupted because 
of the greatly enhanced Rayleigh- Taylor instability. The fragments axe likely to form 
stars. Radio galaxies bigger than 50 kpc are suggested to be in that phase. In every 
phase, a cooling flow may exist that may form an infalling dense shell immediately 
outside of the wind's bow shock. 
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and also possibly on the red wing due to an infalling cooling flow. Eventually, 
the jet plasma will fill all of the cavity (Fig. Id). The enhanced pressure will 
accelerate the bow shock of the wind, which merges with the jet's bow shock 
at that time. This causes the shell to be highly Rayleigh- Taylor unstable, and 
the shell will be disrupted quickly. Additionally, the beams pierce it. The shell 
fragments will fly around in a turbulent way and are likely to form stars. 

This scenario was examined by hydrodynamic simulations, which is re- 
ported on in the following. 



3 Computational aspects 

3.1 Numerics 

For the computations in this contribution, the magneto-hydrodynamic (MHD) 
code Nirvana was employed (Ziegler & Yorke 1997). In that version, it solves 
the MHD equations in three dimensions (3D) for density p, velocity v, internal 
energy e, and magnetic field B: 

(1) 
-Vp-pV^> + -J-(B- V)B--J-VB2 (2) 

47r OTT 

-pV-v-p'^A (3) 

V X (v X B) , (4) 

where ^ denotes an external gravitational potential and yl is a temperature 
dependent cooling fimction (zero metals, zero photoionising field) according to 
Sutherland & Dopita (1993). NIRVANA can be characterised by the following 
properties: 

1. explicit Eulerian time-stepping, 

2. operator-splitting formalism for the advection part of the solver, 

3. method of characteristics-constraint-transport algorithm to solve the in- 
duction equation and to compute the Lorentz forces; 

4. artificial viscosity has been included to dissipate high-frequency noise and 

to allow for shock smearing in case the fiow becomes supersonic. 

The code was vectorised and parallelised by OpenMP like methods, and 
successfully run on the SX-5 (Krause & Camenzind 2002). All the significant 
loops could be vectorised. The number crunching part scales without signifi- 
cant performance loss. This is also true for the MHD part of the solver. The 
code was now modified slightly in order to run on the new NEC SX-6 installa- 
tion of the HLRS. A profile output for a test run is shown in table 1. The test 
run evolved a hydrodynamic problem on a 4096 x 4000 grid for 100 timesteps. 
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Table 1. Profile output: SX-6 test run 





1Q9 SQ7Qfl7 


TT^PT T^lTTIP ( QPr* 1 


771 848942 


Slvt; T'lmp (tjppl 




Vpptnr T'lmp (tjopi 


672 636268 


Xlloli ■ V-vL* UllL 




V Trmt rJoiint 

V . Xilol). \_y>JLilllj 


2366061 1 3nn 


\l P.lpmpnf r^r^nnf 




FT.OP Crmnt 


240401 Q2^748Q 


MOPS 


7888 fi,'in237 


MFLOPS 


3231 227151 


MOP*-i ('rnnriirrpnf ^ 


^71 1 Qf) 




2341 1 01 3720 


A V T.ono'tTi 

iT. . V . J-J'^iltlLli 


2^^ 84032"! 


V On Ratio (%\ 


99 420335 


Memory Size (MB) 


2320.000000 


Mcix Concurrent Proc. 


8 


Cone. Time(>= 1) (sec) 


106.531878 


Cone. Time(>= 2) (sec) 


95.892731 


Cone. Time(>= 3) (sec) 


95.830364 


Cone. Time(>= 4) (see) 


95.813671 


Cone. Time(>= 5) (sec) 


95.798599 


Cone. Timc(>= 6) (sec) 


95.782837 


Cone. Time(>= 7) (sec) 


95.603244 


Cone. Timc(>= 8) (sec) 


90.646246 


Event Busy Count 





Event Wait (sec) 


0.000000 


Lock Busy Count 


17078 


Lock Wait (sec) 


5.299937 


Barrier Busy Count 





Barrier Wait (see) 


0.000000 


MIPS 


76.382190 


MIPS (concurrent) 


553.407238 


I-Cache (see) 


0.408468 


O-Cache (sec) 


6.734071 


Bank (sec) 


21.394816 







This output indicates a performance of roughly 23 GFfop for a run with 8 
processors and an acceptable load ballancing. A single processor run yielded 
3.5 GFlops. This is 38% of the peak performance of 9.2 GFlops. A critical 
point in achieving this quite good performance was a change in the upwind 
part of the advection solver. Now, both directions are computed, and the one 
needed is selected afterwards. This saves an if-clause in the vectorised loops. 
Although, more floating point operations are performed, the speed increase is 
sufficiently high. The code scales well with the number of processors within 
one node (Fig. 2). 

We report to computations. An axisymmetric simulation has been per- 
formed on one processor of the SX-5, involving about 200 CPU hours and 
200 MB of memory. Also, a fully 3D computation has been carried out on two 
processors of the SX-6, involving 100 CPU hours and 2 GB of memory. The 
main reason for the short computation time is the limit of 100 CPU hours for 
the batch queue on the SX-6 and the fact that the machine was installed only 
two weeks ago at the time of writing. We typically need some 1000 CPU hours 
for a full run. Also this run is currently repeated on the SX-5 and evolved for 
a longer time. 

3.2 Setup 

We start the computations with an isothermal (10^ K) King profile for the 
initial density p^: 
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Fig. 2. Scaling of the NIRVANA code on tlie new SX-6 installation at HLRS 



0.3mr 



lOkpc 



(5) 



where r is the spherical radius. The atmosphere is stabilised by the background 
gravity of the dark matter halo in hydrostatic equilibrium. We start the galac- 
tic wind by injecting mass at a rate M — IOMq/jt and thermal energy at a 
rate Lg = 10^^ erg/yr, which corresponds to one supernova per year, inside a 
region of 3 kpc, with an exponential decay. For the 2.5D run (axisymmetry), 
the grid has 2047 x 1023 points that correspond to 200 x 100 kpc. In the 3D 
case, the grid has 511 x 201 x 201 points, corresponding to 200 x 80 x 80 kpc. 
With a jet radius of one (two) kpc the resolution corresponds to 10 (5) points 
per beam radius for the 2.5D (3D) case. In both cases, the resolution is enough 
to resolve the radiative bow shock. 

After 80 Myrs, when the dense shell has formed, bipolar jets are injected 
in the center of the galaxy which also corresponds to the center of the grid. 
The jet density is 10~^mp cm~^. This is roughly a hundred times less than 
the gas inside of the wind bubble at the time of the jet start. The jets are 
then evolved for 18 (2) Myrs. 
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4.1 Axisymmetric simulation 
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Density and temperature maps of the 2.5D simulation are shown at two 
different times in Fig 3. In the earher plots, the jet is well inside the spherically 
symmetric galactic wind. This wind has already formed a dense shell before 
the start of the jet. Careful inspection of Fig 3c reveals the shell to have a 
temperature of 10^ K, with a thin enhancement to « 10^ K in the middle. 
This is the position of the bow shock. The structure can be seen in more detail 
in Fig 4. The bow shock of the galactic wind expands into a cooling flow. The 
shock, at 400 km/s, heats the gas to > 10^ K. Inside of the shock, the gas cools 
quickly down to 10* K, the lower end of the cooling function. Outside of the 
shock, the cooling flow has also produced a cool and dense shell, inflowing at 
« 100 km/s. The two jets pierce a considerable hole into these shells. When the 
jet's bow shock hits the one of the wind, the shell gets accelerated. This shell 
can only remain stable as long as its deceleration exceeds the local gravity. 
An acceleration greatly enhances the Rayleigh- Taylor instability. Hence, a 
considerable amount of gas is entrained. Fig. 3b shows this entrained gas in 
yellow-green colors, extending in fingers from the shell down to the jet beam, 
where it assembles due to gravity. The shell is transformed into clumps that 
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Fig. 4. Slice of density, temperature and sideways velocity for the 2.5D run at 
80.5 Myr. The slice is located vertical to the jet in the center. The bow shock is 
located at = 29 kpc. Inside of it the collapsed shell flows outward at « 400 km/s, 
outside of it cool gas flows in at fs 100 km/s. 

still get denser at the end of the simulation. A magnification of the bow 
shock region in Fig. 3d is shown in Fig. 5. The fragments of the shell can be 
seen in blue. Between the fragments, the very hot (red) jet plasma, that was 
injected through the beams and distributed in the wind cavity, flows past and 
compresses the ambient gas, which drives a faster, non-radiative bow shock. 

The dissolution of the shell can also be seen from the number density 
versus sideways velocity histograms (Fig 6). At 80.5 Myr, the motion of the 
dense gas is well ordered, centered at 200 km/s. At 95 Myr, higher densities 
and higher velocities are present, but the distribution is much less uniform. 
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Fig. 5. Magnification of the bow shock region of Fig. 3d. The highest temperature 
has the exhaust plasma from the jet beam. This plasma (hot, red) flows through the 
holes between the fragments (cold, blue) of the shell and compresses the ambient 
gas thereby establishing a new, non-radiative bow shock. The highest temperature 
is shown in white. 




Fig. 6. Frequency of number densities and sideways (radial) velocities at 80.5 Myr 
and 95 Myr. The counts are cut at 10^. 
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4.2 Three-dimensional simulation 




Fig. 7. Density (top) and Mach number in jet direction (bottom) for the 3D simu- 
lation at 82 Myr. 



The jets in the 3D run stay quite symmetric. Only the cocoon shows 
some three-dimensional structure. The bow shock of the jet has unfortunately 
not yet reached the wind's bow shock, and the process of fragmentation can 
therefore not yet be compared to the 2.5D case. An interesting result is that 
the Mach number drops considerably at the location of the shell. 



5 Discussion 



Interaction of Jets with Galactic Winds 



11 



The simulations show that a jet is able to disrupt a wind blown shell via 
the Rayleigh- Taylor instability during the simulation time. This was shown 
directly in density contour plots, as well as the histograms of velocity versus 
number density. The 3D computation has not yet been evolved for long enough 
in order to compare the fragmentation process. Another interesting feature to 
appear in 3D is jet bending. Slight changes of the jet direction might cause the 
beam to hit the shell on one side, which may deflect it. This effect is absent, the 
beams stay straight. An emission line region within such a bubble would first 
be absorbed on the blue wing, and if a cooling flow is present, possibly also 
on the red wing. Observed radio galaxies in the young universe are absorbed 
preferentially on the blue wing. This indicates that, if the suggested scenario is 
correct, the gas around most of these objects is cither less dense or has a higher 
temperature, in order not to show a cooling flow. However, the properties of 
the environments of these objects seem to be quite uniform (van Ojik et al. 
1997), because of the strict absence of absorbers for objects larger than 50 kpc, 
and the presence in nearly all of the objects smaller than that. It is therefore 
likely that they are located at the border line. An inflowing cool shell will be 
created when the cooling time for the ambient gas is shorter than the delay 
between starburst and jet start (« 10® years): 



For the gas inside of the wind's bow shock to cool, the shock's Macli number 
should be close to unity. Since absorbers are typically blue-shifted by a few 
hundred km/s, it follows: 



These two requirements fit together, if the ambient densities arc roughly 
0.1 cm~^, and the temperature is typically about 10^ K, i.e. the galaxy clus- 
ters in the young universe would have had denser and colder gas than nearby 
ones. 

This comparison, and the simulation results achieved on the NEC SX- 
5 and SX 6 supercomputers at the HLRS suggest that high redshift radio 
galaxies may indeed be associated with thin shells blown by a galactic wind. 
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10^ yr = tcooi.ambient = 12 Myr VT/107K/(n/cm-3) 



(6) 



200 km/s « Csound.ambient = 400 km/s s/tJWk 



(7) 
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